home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / quadpack / dqagpe.f < prev    next >
Text File  |  1996-07-19  |  21KB  |  557 lines
  1.       SUBROUTINE DQAGPE(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,LIMIT,RESULT,
  2.      *   ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,PTS,IORD,LEVEL,NDIN,
  3.      *   LAST)
  4. C***BEGIN PROLOGUE  DQAGPE
  5. C***DATE WRITTEN   800101   (YYMMDD)
  6. C***REVISION DATE  830518   (YYMMDD)
  7. C***CATEGORY NO.  H2A2A1
  8. C***KEYWORDS  AUTOMATIC INTEGRATOR, GENERAL-PURPOSE,
  9. C             SINGULARITIES AT USER SPECIFIED POINTS,
  10. C             EXTRAPOLATION, GLOBALLY ADAPTIVE.
  11. C***AUTHOR  PIESSENS,ROBERT ,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  12. C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  13. C***PURPOSE  THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
  14. C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B), HOPEFULLY
  15. C            SATISFYING FOLLOWING CLAIM FOR ACCURACY ABS(I-RESULT).LE.
  16. C            MAX(EPSABS,EPSREL*ABS(I)). BREAK POINTS OF THE INTEGRATION
  17. C            INTERVAL, WHERE LOCAL DIFFICULTIES OF THE INTEGRAND MAY
  18. C            OCCUR(E.G. SINGULARITIES,DISCONTINUITIES),PROVIDED BY USER.
  19. C***DESCRIPTION
  20. C
  21. C        COMPUTATION OF A DEFINITE INTEGRAL
  22. C        STANDARD FORTRAN SUBROUTINE
  23. C        DOUBLE PRECISION VERSION
  24. C
  25. C        PARAMETERS
  26. C         ON ENTRY
  27. C            F      - DOUBLE PRECISION
  28. C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  29. C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  30. C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  31. C
  32. C            A      - DOUBLE PRECISION
  33. C                     LOWER LIMIT OF INTEGRATION
  34. C
  35. C            B      - DOUBLE PRECISION
  36. C                     UPPER LIMIT OF INTEGRATION
  37. C
  38. C            NPTS2  - INTEGER
  39. C                     NUMBER EQUAL TO TWO MORE THAN THE NUMBER OF
  40. C                     USER-SUPPLIED BREAK POINTS WITHIN THE INTEGRATION
  41. C                     RANGE, NPTS2.GE.2.
  42. C                     IF NPTS2.LT.2, THE ROUTINE WILL END WITH IER = 6.
  43. C
  44. C            POINTS - DOUBLE PRECISION
  45. C                     VECTOR OF DIMENSION NPTS2, THE FIRST (NPTS2-2)
  46. C                     ELEMENTS OF WHICH ARE THE USER PROVIDED BREAK
  47. C                     POINTS. IF THESE POINTS DO NOT CONSTITUTE AN
  48. C                     ASCENDING SEQUENCE THERE WILL BE AN AUTOMATIC
  49. C                     SORTING.
  50. C
  51. C            EPSABS - DOUBLE PRECISION
  52. C                     ABSOLUTE ACCURACY REQUESTED
  53. C            EPSREL - DOUBLE PRECISION
  54. C                     RELATIVE ACCURACY REQUESTED
  55. C                     IF  EPSABS.LE.0
  56. C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
  57. C                     THE ROUTINE WILL END WITH IER = 6.
  58. C
  59. C            LIMIT  - INTEGER
  60. C                     GIVES AN UPPER BOUND ON THE NUMBER OF SUBINTERVALS
  61. C                     IN THE PARTITION OF (A,B), LIMIT.GE.NPTS2
  62. C                     IF LIMIT.LT.NPTS2, THE ROUTINE WILL END WITH
  63. C                     IER = 6.
  64. C
  65. C         ON RETURN
  66. C            RESULT - DOUBLE PRECISION
  67. C                     APPROXIMATION TO THE INTEGRAL
  68. C
  69. C            ABSERR - DOUBLE PRECISION
  70. C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  71. C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
  72. C
  73. C            NEVAL  - INTEGER
  74. C                     NUMBER OF INTEGRAND EVALUATIONS
  75. C
  76. C            IER    - INTEGER
  77. C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
  78. C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
  79. C                             ACCURACY HAS BEEN ACHIEVED.
  80. C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE.
  81. C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
  82. C                             LESS RELIABLE. IT IS ASSUMED THAT THE
  83. C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
  84. C                      IER.LT.0 EXIT REQUESTED FROM USER-SUPPLIED
  85. C                             FUNCTION.
  86. C
  87. C            ERROR MESSAGES
  88. C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
  89. C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
  90. C                             SUBDIVISIONS BY INCREASING THE VALUE OF
  91. C                             LIMIT (AND TAKING THE ACCORDING DIMENSION
  92. C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
  93. C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
  94. C                             TO ANALYZE THE INTEGRAND IN ORDER TO
  95. C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
  96. C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
  97. C                             DETERMINED (I.E. SINGULARITY,
  98. C                             DISCONTINUITY WITHIN THE INTERVAL), IT
  99. C                             SHOULD BE SUPPLIED TO THE ROUTINE AS AN
  100. C                             ELEMENT OF THE VECTOR POINTS. IF NECESSARY
  101. C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
  102. C                             MUST BE USED, WHICH IS DESIGNED FOR
  103. C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
  104. C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
  105. C                             DETECTED, WHICH PREVENTS THE REQUESTED
  106. C                             TOLERANCE FROM BEING ACHIEVED.
  107. C                             THE ERROR MAY BE UNDER-ESTIMATED.
  108. C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
  109. C                             AT SOME POINTS OF THE INTEGRATION
  110. C                             INTERVAL.
  111. C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
  112. C                             ROUNDOFF ERROR IS DETECTED IN THE
  113. C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
  114. C                             THE REQUESTED TOLERANCE CANNOT BE
  115. C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
  116. C                             THE BEST WHICH CAN BE OBTAINED.
  117. C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
  118. C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
  119. C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
  120. C                             OF IER.GT.0.
  121. C                         = 6 THE INPUT IS INVALID BECAUSE
  122. C                             NPTS2.LT.2 OR
  123. C                             BREAK POINTS ARE SPECIFIED OUTSIDE
  124. C                             THE INTEGRATION RANGE OR
  125. C                             (EPSABS.LE.0 AND
  126. C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28))
  127. C                             OR LIMIT.LT.NPTS2.
  128. C                             RESULT, ABSERR, NEVAL, LAST, RLIST(1),
  129. C                             AND ELIST(1) ARE SET TO ZERO. ALIST(1) AND
  130. C                             BLIST(1) ARE SET TO A AND B RESPECTIVELY.
  131. C
  132. C            ALIST  - DOUBLE PRECISION
  133. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  134. C                      LAST  ELEMENTS OF WHICH ARE THE LEFT END POINTS
  135. C                     OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
  136. C                     INTEGRATION RANGE (A,B)
  137. C
  138. C            BLIST  - DOUBLE PRECISION
  139. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  140. C                      LAST  ELEMENTS OF WHICH ARE THE RIGHT END POINTS
  141. C                     OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
  142. C                     INTEGRATION RANGE (A,B)
  143. C
  144. C            RLIST  - DOUBLE PRECISION
  145. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  146. C                      LAST  ELEMENTS OF WHICH ARE THE INTEGRAL
  147. C                     APPROXIMATIONS ON THE SUBINTERVALS
  148. C
  149. C            ELIST  - DOUBLE PRECISION
  150. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
  151. C                      LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
  152. C                     ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
  153. C
  154. C            PTS    - DOUBLE PRECISION
  155. C                     VECTOR OF DIMENSION AT LEAST NPTS2, CONTAINING THE
  156. C                     INTEGRATION LIMITS AND THE BREAK POINTS OF THE
  157. C                     INTERVAL IN ASCENDING SEQUENCE.
  158. C
  159. C            LEVEL  - INTEGER
  160. C                     VECTOR OF DIMENSION AT LEAST LIMIT, CONTAINING THE
  161. C                     SUBDIVISION LEVELS OF THE SUBINTERVAL, I.E. IF
  162. C                     (AA,BB) IS A SUBINTERVAL OF (P1,P2) WHERE P1 AS
  163. C                     WELL AS P2 IS A USER-PROVIDED BREAK POINT OR
  164. C                     INTEGRATION LIMIT, THEN (AA,BB) HAS LEVEL L IF
  165. C                     ABS(BB-AA) = ABS(P2-P1)*2**(-L).
  166. C
  167. C            NDIN   - INTEGER
  168. C                     VECTOR OF DIMENSION AT LEAST NPTS2, AFTER FIRST
  169. C                     INTEGRATION OVER THE INTERVALS (PTS(I)),PTS(I+1),
  170. C                     I = 0,1, ..., NPTS2-2, THE ERROR ESTIMATES OVER
  171. C                     SOME OF THE INTERVALS MAY HAVE BEEN INCREASED
  172. C                     ARTIFICIALLY, IN ORDER TO PUT THEIR SUBDIVISION
  173. C                     FORWARD. IF THIS HAPPENS FOR THE SUBINTERVAL
  174. C                     NUMBERED K, NDIN(K) IS PUT TO 1, OTHERWISE
  175. C                     NDIN(K) = 0.
  176. C
  177. C            IORD   - INTEGER
  178. C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
  179. C                     ELEMENTS OF WHICH ARE POINTERS TO THE
  180. C                     ERROR ESTIMATES OVER THE SUBINTERVALS,
  181. C                     SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
  182. C                     FORM A DECREASING SEQUENCE, WITH K = LAST
  183. C                     IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
  184. C                     OTHERWISE
  185. C
  186. C            LAST   - INTEGER
  187. C                     NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
  188. C                     SUBDIVISIONS PROCESS
  189. C
  190. C***REFERENCES  (NONE)
  191. C***ROUTINES CALLED  D1MACH,DQELG,DQK21,DQPSRT
  192. C***END PROLOGUE  DQAGPE
  193.       DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
  194.      *  A2,B,BLIST,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2,DMAX1,DMIN1,
  195.      *  DRES,D1MACH,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,
  196.      *  ERRMAX,ERROR1,ERRO12,ERROR2,ERRSUM,ERTEST,F,OFLOW,POINTS,PTS,
  197.      *  RESA,RESABS,RESEPS,RESULT,RES3LA,RLIST,RLIST2,SIGN,TEMP,UFLOW
  198.       INTEGER I,ID,IER,IERRO,IND1,IND2,IORD,IP1,IROFF1,IROFF2,IROFF3,J,
  199.      *  JLOW,JUPBND,K,KSGN,KTMIN,LAST,LEVCUR,LEVEL,LEVMAX,LIMIT,MAXERR,
  200.      *  NDIN,NEVAL,NINT,NINTP1,NPTS,NPTS2,NRES,NRMAX,NUMRL2
  201.       LOGICAL EXTRAP,NOEXT
  202. C
  203. C
  204.       DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
  205.      *  LEVEL(LIMIT),NDIN(NPTS2),POINTS(NPTS2),PTS(NPTS2),RES3LA(3),
  206.      *  RLIST(LIMIT),RLIST2(52)
  207. C
  208.       EXTERNAL F
  209. C
  210. C            THE DIMENSION OF RLIST2 IS DETERMINED BY THE VALUE OF
  211. C            LIMEXP IN SUBROUTINE EPSALG (RLIST2 SHOULD BE OF DIMENSION
  212. C            (LIMEXP+2) AT LEAST).
  213. C
  214. C
  215. C            LIST OF MAJOR VARIABLES
  216. C            -----------------------
  217. C
  218. C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
  219. C                       CONSIDERED UP TO NOW
  220. C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
  221. C                       CONSIDERED UP TO NOW
  222. C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
  223. C                       (ALIST(I),BLIST(I))
  224. C           RLIST2    - ARRAY OF DIMENSION AT LEAST LIMEXP+2
  225. C                       CONTAINING THE PART OF THE EPSILON TABLE WHICH
  226. C                       IS STILL NEEDED FOR FURTHER COMPUTATIONS
  227. C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
  228. C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST ERROR
  229. C                       ESTIMATE
  230. C           ERRMAX    - ELIST(MAXERR)
  231. C           ERLAST    - ERROR ON THE INTERVAL CURRENTLY SUBDIVIDED
  232. C                       (BEFORE THAT SUBDIVISION HAS TAKEN PLACE)
  233. C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
  234. C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
  235. C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
  236. C                       ABS(RESULT))
  237. C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
  238. C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
  239. C           LAST      - INDEX FOR SUBDIVISION
  240. C           NRES      - NUMBER OF CALLS TO THE EXTRAPOLATION ROUTINE
  241. C           NUMRL2    - NUMBER OF ELEMENTS IN RLIST2. IF AN APPROPRIATE
  242. C                       APPROXIMATION TO THE COMPOUNDED INTEGRAL HAS
  243. C                       BEEN OBTAINED, IT IS PUT IN RLIST2(NUMRL2) AFTER
  244. C                       NUMRL2 HAS BEEN INCREASED BY ONE.
  245. C           ERLARG    - SUM OF THE ERRORS OVER THE INTERVALS LARGER
  246. C                       THAN THE SMALLEST INTERVAL CONSIDERED UP TO NOW
  247. C           EXTRAP    - LOGICAL VARIABLE DENOTING THAT THE ROUTINE
  248. C                       IS ATTEMPTING TO PERFORM EXTRAPOLATION. I.E.
  249. C                       BEFORE SUBDIVIDING THE SMALLEST INTERVAL WE
  250. C                       TRY TO DECREASE THE VALUE OF ERLARG.
  251. C           NOEXT     - LOGICAL VARIABLE DENOTING THAT EXTRAPOLATION IS
  252. C                       NO LONGER ALLOWED (TRUE-VALUE)
  253. C
  254. C            MACHINE DEPENDENT CONSTANTS
  255. C            ---------------------------
  256. C
  257. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  258. C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  259. C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
  260. C
  261. C***FIRST EXECUTABLE STATEMENT  DQAGPE
  262.       EPMACH = D1MACH(4)
  263. C
  264. C            TEST ON VALIDITY OF PARAMETERS
  265. C            -----------------------------
  266. C
  267.       IER = 0
  268.       NEVAL = 0
  269.       LAST = 0
  270.       RESULT = 0.0D+00
  271.       ABSERR = 0.0D+00
  272.       ALIST(1) = A
  273.       BLIST(1) = B
  274.       RLIST(1) = 0.0D+00
  275.       ELIST(1) = 0.0D+00
  276.       IORD(1) = 0
  277.       LEVEL(1) = 0
  278.       NPTS = NPTS2-2
  279.       IF(NPTS2.LT.2.OR.LIMIT.LE.NPTS.OR.(EPSABS.LE.0.0D+00.AND.
  280.      *  EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28))) IER = 6
  281.       IF(IER.EQ.6) GO TO 999
  282. C
  283. C            IF ANY BREAK POINTS ARE PROVIDED, SORT THEM INTO AN
  284. C            ASCENDING SEQUENCE.
  285. C
  286.       SIGN = 1.0D+00
  287.       IF(A.GT.B) SIGN = -1.0D+00
  288.       PTS(1) = DMIN1(A,B)
  289.       IF(NPTS.EQ.0) GO TO 15
  290.       DO 10 I = 1,NPTS
  291.         PTS(I+1) = POINTS(I)
  292.    10 CONTINUE
  293.    15 PTS(NPTS+2) = DMAX1(A,B)
  294.       NINT = NPTS+1
  295.       A1 = PTS(1)
  296.       IF(NPTS.EQ.0) GO TO 40
  297.       NINTP1 = NINT+1
  298.       DO 20 I = 1,NINT
  299.         IP1 = I+1
  300.         DO 20 J = IP1,NINTP1
  301.           IF(PTS(I).LE.PTS(J)) GO TO 20
  302.           TEMP = PTS(I)
  303.           PTS(I) = PTS(J)
  304.           PTS(J) = TEMP
  305.    20 CONTINUE
  306.       IF(PTS(1).NE.DMIN1(A,B).OR.PTS(NINTP1).NE.DMAX1(A,B)) IER = 6
  307.       IF(IER.EQ.6) GO TO 999
  308. C
  309. C            COMPUTE FIRST INTEGRAL AND ERROR APPROXIMATIONS.
  310. C            ------------------------------------------------
  311. C
  312.    40 RESABS = 0.0D+00
  313.       DO 50 I = 1,NINT
  314.         B1 = PTS(I+1)
  315.         CALL DQK21(F,A1,B1,AREA1,ERROR1,DEFABS,RESA,IER)
  316.         IF (IER .LT. 0) RETURN
  317.         ABSERR = ABSERR+ERROR1
  318.         RESULT = RESULT+AREA1
  319.         NDIN(I) = 0
  320.         IF(ERROR1.EQ.RESA.AND.ERROR1.NE.0.0D+00) NDIN(I) = 1
  321.         RESABS = RESABS+DEFABS
  322.         LEVEL(I) = 0
  323.         ELIST(I) = ERROR1
  324.         ALIST(I) = A1
  325.         BLIST(I) = B1
  326.         RLIST(I) = AREA1
  327.         IORD(I) = I
  328.         A1 = B1
  329.    50 CONTINUE
  330.       ERRSUM = 0.0D+00
  331.       DO 55 I = 1,NINT
  332.         IF(NDIN(I).EQ.1) ELIST(I) = ABSERR
  333.         ERRSUM = ERRSUM+ELIST(I)
  334.    55 CONTINUE
  335. C
  336. C           TEST ON ACCURACY.
  337. C
  338.       LAST = NINT
  339.       NEVAL = 21*NINT
  340.       DRES = DABS(RESULT)
  341.       ERRBND = DMAX1(EPSABS,EPSREL*DRES)
  342.       IF(ABSERR.LE.0.1D+03*EPMACH*RESABS.AND.ABSERR.GT.ERRBND) IER = 2
  343.       IF(NINT.EQ.1) GO TO 80
  344.       DO 70 I = 1,NPTS
  345.         JLOW = I+1
  346.         IND1 = IORD(I)
  347.         DO 60 J = JLOW,NINT
  348.           IND2 = IORD(J)
  349.           IF(ELIST(IND1).GT.ELIST(IND2)) GO TO 60
  350.           IND1 = IND2
  351.           K = J
  352.    60   CONTINUE
  353.         IF(IND1.EQ.IORD(I)) GO TO 70
  354.         IORD(K) = IORD(I)
  355.         IORD(I) = IND1
  356.    70 CONTINUE
  357.       IF(LIMIT.LT.NPTS2) IER = 1
  358.    80 IF(IER.NE.0.OR.ABSERR.LE.ERRBND) GO TO 210
  359. C
  360. C           INITIALIZATION
  361. C           --------------
  362. C
  363.       RLIST2(1) = RESULT
  364.       MAXERR = IORD(1)
  365.       ERRMAX = ELIST(MAXERR)
  366.       AREA = RESULT
  367.       NRMAX = 1
  368.       NRES = 0
  369.       NUMRL2 = 1
  370.       KTMIN = 0
  371.       EXTRAP = .FALSE.
  372.       NOEXT = .FALSE.
  373.       ERLARG = ERRSUM
  374.       ERTEST = ERRBND
  375.       LEVMAX = 1
  376.       IROFF1 = 0
  377.       IROFF2 = 0
  378.       IROFF3 = 0
  379.       IERRO = 0
  380.       UFLOW = D1MACH(1)
  381.       OFLOW = D1MACH(2)
  382.       ABSERR = OFLOW
  383.       KSGN = -1
  384.       IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*RESABS) KSGN = 1
  385. C
  386. C           MAIN DO-LOOP
  387. C           ------------
  388. C
  389.       DO 160 LAST = NPTS2,LIMIT
  390. C
  391. C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
  392. C           ESTIMATE.
  393. C
  394.         LEVCUR = LEVEL(MAXERR)+1
  395.         A1 = ALIST(MAXERR)
  396.         B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
  397.         A2 = B1
  398.         B2 = BLIST(MAXERR)
  399.         ERLAST = ERRMAX
  400.         CALL DQK21(F,A1,B1,AREA1,ERROR1,RESA,DEFAB1,IER)
  401.         IF (IER .LT. 0) RETURN
  402.         CALL DQK21(F,A2,B2,AREA2,ERROR2,RESA,DEFAB2,IER)
  403.         IF (IER .LT. 0) RETURN
  404. C
  405. C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
  406. C           AND ERROR AND TEST FOR ACCURACY.
  407. C
  408.         NEVAL = NEVAL+42
  409.         AREA12 = AREA1+AREA2
  410.         ERRO12 = ERROR1+ERROR2
  411.         ERRSUM = ERRSUM+ERRO12-ERRMAX
  412.         AREA = AREA+AREA12-RLIST(MAXERR)
  413.         IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 95
  414.         IF(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12)
  415.      *  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 90
  416.         IF(EXTRAP) IROFF2 = IROFF2+1
  417.         IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
  418.    90   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
  419.    95   LEVEL(MAXERR) = LEVCUR
  420.         LEVEL(LAST) = LEVCUR
  421.         RLIST(MAXERR) = AREA1
  422.         RLIST(LAST) = AREA2
  423.         ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA))
  424. C
  425. C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
  426. C
  427.         IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
  428.         IF(IROFF2.GE.5) IERRO = 3
  429. C
  430. C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF
  431. C           SUBINTERVALS EQUALS LIMIT.
  432. C
  433.         IF(LAST.EQ.LIMIT) IER = 1
  434. C
  435. C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
  436. C           AT A POINT OF THE INTEGRATION RANGE
  437. C
  438.         IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
  439.      *  (DABS(A2)+0.1D+04*UFLOW)) IER = 4
  440. C
  441. C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
  442. C
  443.         IF(ERROR2.GT.ERROR1) GO TO 100
  444.         ALIST(LAST) = A2
  445.         BLIST(MAXERR) = B1
  446.         BLIST(LAST) = B2
  447.         ELIST(MAXERR) = ERROR1
  448.         ELIST(LAST) = ERROR2
  449.         GO TO 110
  450.   100   ALIST(MAXERR) = A2
  451.         ALIST(LAST) = A1
  452.         BLIST(LAST) = B1
  453.         RLIST(MAXERR) = AREA2
  454.         RLIST(LAST) = AREA1
  455.         ELIST(MAXERR) = ERROR2
  456.         ELIST(LAST) = ERROR1
  457. C
  458. C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
  459. C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
  460. C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
  461. C
  462.   110   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
  463. C ***JUMP OUT OF DO-LOOP
  464.         IF(ERRSUM.LE.ERRBND) GO TO 190
  465. C ***JUMP OUT OF DO-LOOP
  466.         IF(IER.NE.0) GO TO 170
  467.         IF(NOEXT) GO TO 160
  468.         ERLARG = ERLARG-ERLAST
  469.         IF(LEVCUR+1.LE.LEVMAX) ERLARG = ERLARG+ERRO12
  470.         IF(EXTRAP) GO TO 120
  471. C
  472. C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
  473. C           SMALLEST INTERVAL.
  474. C
  475.         IF(LEVEL(MAXERR)+1.LE.LEVMAX) GO TO 160
  476.         EXTRAP = .TRUE.
  477.         NRMAX = 2
  478.   120   IF(IERRO.EQ.3.OR.ERLARG.LE.ERTEST) GO TO 140
  479. C
  480. C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
  481. C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER
  482. C           THE LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
  483. C
  484.         ID = NRMAX
  485.         JUPBND = LAST
  486.         IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
  487.         DO 130 K = ID,JUPBND
  488.           MAXERR = IORD(NRMAX)
  489.           ERRMAX = ELIST(MAXERR)
  490. C ***JUMP OUT OF DO-LOOP
  491.           IF(LEVEL(MAXERR)+1.LE.LEVMAX) GO TO 160
  492.           NRMAX = NRMAX+1
  493.   130   CONTINUE
  494. C
  495. C           PERFORM EXTRAPOLATION.
  496. C
  497.   140   NUMRL2 = NUMRL2+1
  498.         RLIST2(NUMRL2) = AREA
  499.         IF(NUMRL2.LE.2) GO TO 155
  500.         CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
  501.         KTMIN = KTMIN+1
  502.         IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
  503.         IF(ABSEPS.GE.ABSERR) GO TO 150
  504.         KTMIN = 0
  505.         ABSERR = ABSEPS
  506.         RESULT = RESEPS
  507.         CORREC = ERLARG
  508.         ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS))
  509. C ***JUMP OUT OF DO-LOOP
  510.         IF(ABSERR.LT.ERTEST) GO TO 170
  511. C
  512. C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
  513. C
  514.   150   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
  515.         IF(IER.GE.5) GO TO 170
  516.   155   MAXERR = IORD(1)
  517.         ERRMAX = ELIST(MAXERR)
  518.         NRMAX = 1
  519.         EXTRAP = .FALSE.
  520.         LEVMAX = LEVMAX+1
  521.         ERLARG = ERRSUM
  522.   160 CONTINUE
  523. C
  524. C           SET THE FINAL RESULT.
  525. C           ---------------------
  526. C
  527. C
  528.   170 IF(ABSERR.EQ.OFLOW) GO TO 190
  529.       IF((IER+IERRO).EQ.0) GO TO 180
  530.       IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
  531.       IF(IER.EQ.0) IER = 3
  532.       IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00)GO TO 175
  533.       IF(ABSERR.GT.ERRSUM)GO TO 190
  534.       IF(AREA.EQ.0.0D+00) GO TO 210
  535.       GO TO 180
  536.   175 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA))GO TO 190
  537. C
  538. C           TEST ON DIVERGENCE.
  539. C
  540.   180 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE.
  541.      *  RESABS*0.1D-01) GO TO 210
  542.       IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03.OR.
  543.      *  ERRSUM.GT.DABS(AREA)) IER = 6
  544.       GO TO 210
  545. C
  546. C           COMPUTE GLOBAL INTEGRAL SUM.
  547. C
  548.   190 RESULT = 0.0D+00
  549.       DO 200 K = 1,LAST
  550.         RESULT = RESULT+RLIST(K)
  551.   200 CONTINUE
  552.       ABSERR = ERRSUM
  553.   210 IF(IER.GT.2) IER = IER-1
  554.       RESULT = RESULT*SIGN
  555.   999 RETURN
  556.       END
  557.